Lattice Boltzmann Simulations of Liquid Crystal Hydrodynamics 
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We describe a lattice Boltzmann algorithm to simulate liquid crystal hydrodynamics. The equa- 
tions of motion are written in terms of a tensor order parameter. This allows both the isotropic and 
the nematic phases to be considered. Backflow effects and the hydrodynamics of topological defects 
are naturally included in the simulations, as are viscoelastic properties such as shear-thinning and 
shear-banding. 
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I. INTRODUCTION 



Liquid crystalline materials are often made up of long, thin, rod-like molecules The molecular geometry and 
^ . interactions can lead to a wide range of equilibrium phases. Here we shall be concerned with two of the simplest, the 
isotropic phase, where the orientation of the molecules is random, and the nematic phase, where the molecules tend 
to align along a preferred direction. 

The aim of this paper is to describe a numerical scheme which can explore the hydrodynamics of liquid crystals 
| within both the isotropic and the nematic phases. There are two major differences between the hydrodynamics of 
<—j . simple liquids and that of liquid crystals. First, the geometry of the molecules means that they are rotated by gradients 
in the velocity field. Second, the equilibrium free energy is more complex than for a simple fluid and this in turn 
increases the complexity of the stress tensor in the Navier-Stokes equation for the evolution of the fluid momentum. 
This coupling between the elastic energy and the flow leads to rich hydrodynamic behaviour. A simple example is 
7-H ! the existence of a tumbling phase where the molecules rotate in an applied shear ^ . Other examples include shear 
banding, a non-equilibrium phase separation into coexisting states with different strain rates Q, and the possibility 
of Williams domains, convection cells induced by an applied electric field 0. 

The equations of motion describing liquid crystal hydrodynamics are complex. There are several derivations broadly 
in agreement, but differing in the detailed form of some terms. Here we follow the approach of Beris and Edwards 
Pj who write the equations of motion in terms of a tensor order parameter Q which can be related to the second 
moment of the orientational distribution function of the molecules. This has the advantage that the hydrodynamics 
of both the isotropic and the nematic phases, and of topological defects in the nematic phase, can be included within 
the same formalism. Most other theories of liquid crystal hydrodynamics appear as limiting cases. In particular the 
Ericksen-Leslie formulation of nematodynamics widely used in the experimental liquid crystal literature, follows 
when uniaxiality is imposed and the magnitude of the order parameter is held constant. 

Considerable analytic progress in understanding liquid crystal flow in simple geometries has been made, but this 
is inevitably limited by the complexity of the equations of motion. Therefore it is useful to formulate a method of 
obtaining numerical solutions of the hydrodynamic equations to further explore their rich phenomenology. Moreover we 
should like to be able to predict flow patterns for given viscous and elastic coefficients for comparison to experiments 
and to explore the effects of hydrodynamics when liquid crystals are used in display devices or during industrial 
processing. 

Rey and Tsuji have obtained interesting results on flow-induced ordering of the director field and on defect 
dynamics by solving the Beris-Edwards equation for the order parameter. However, the velocity field was imposed 
externally and no back-flows (effect of the director configuration on the velocity field) were included. Fukuda j?J used 
an Euler scheme to solve a model somewhat simpler than the full Beris-Edwards model but still including backflow, 
and studied the effect of hydrodynamics on phase ordering in liquid crystals. Otherwise most previous work on liquid 
crystal hydrodynamics has been limited to a constant order parameter (the Ericksen-Leslie-Parodi equations) and 
often restricted to one dimension. 

Lattice Boltzmann schemes have recently proved very successful in simulations of complex fluids and it is this 
approach that we shall take here ||. Such algorithms can be usefully and variously considered as a slightly unusual 
finite-difference discretization of the equations of motion or as a lattice version of a simplified Boltzmann equation. It 
is not understood why the approach is particularly useful for complex fluids but it may be related to the very natural 
way in which a free energy describing the equilibrium properties of the fluid can be incorporated in the simulations, 
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drawing on ideas from statistical mechanics [|9|. Recent applications have included phase ordering and flow in binary 
fluids (lCj and self-assembly and spontaneous emulsification in amphiphilic fluids [fLf]|lJ] . 

However, in applications so far, with the exception of |ll|, the order parameter has been a scalar and has coupled 
to the flow via a simple advcctive term. The liquid crystal equations of motion are written in terms of a tensor order 
parameter. This is responsible for the main new features of the lattice Boltzmann approach described in this paper. It 
also leads to the possibility of exploring viscoelastic fluid behaviour such as shear-thinning and shear-banding without 
the need to impose a constitutive equation for the stress JL4| . 

In Section 2 we summarise the hydrodynamic equations of motion for liquid crystals. The lattice Boltzmann 
scheme is defined in Section 3. A modified version of the collision operator is used to eliminate lattice viscosity effects. 
Section 4 describes a Chapman-Enskog expansion which relates the numerical scheme to the hydrodynamic equations 
of motion. Numerical results for simple shear flows are presented in Section 5 and other possible applications of the 
approach are outlined in Section 6. 



II. THE HYDRODYNAMIC EQUATIONS OF MOTION 

We shall follow the formulation of liquid crystal hydrodynamics described by Beris and Edwards [Q . The continuum 
equations of motion are written in terms of a tensor order parameter Q which is related to the direction of individual 
molecules n by Q a p = (n a np — \8 a p) where the angular brackets denote a coarse-grained average. (Greek indices will 
be used to represent Cartesian components of vectors and tensors and the usual summation over repeated indices will 
be assumed.) Q is a traceless symmetric tensor which is zero in the isotropic phase. We first write down a Landau 
free energy which describes the equilibrium properties of the liquid crystal and the isotropic-nematic transition. This 
appears in the equation of motion of the order parameter, which includes a Cahn-Hilliard-like term through which 
the system evolves towards thermodynamic equilibrium. It also includes a term coupling the order parameter to the 
flow. The order parameter is both advected by the flow and, because liquid crystal molecules are rod-like, rotated by 
velocity gradients. 

We then write down the continuity and Navier-Stokes equations for the evolution of the flow field. In particular 
the form of the stress appropriate to a tensor order parameter is discussed. A brief comparison is given to a similar 
formalism introduced by Doi and extended by Olmsted et. al. j[(|[l7). For a uniaxial nematic in the absence 
of any defects the Beris-Edwards equations reduce to the Ericksen-Leslie-Parodi formulation of nematodynamics [Q . 
The hydrodynamic behaviour of nematic liquid crystals is often characterised in terms of the Leslie coefficients and 
it is therefore useful to list them below. More details of the mapping between the Beris-Edwards and the Ericksen- 
Leslie-Parodi equations are given in Appendix A. 

Free energy: The equilibrium properties of a liquid crystal in solution can be described by a free energy |L7| 



F= / d 3 r{ \Ql f} - h -Q a pQp 1 Q ia + ^(Q^f + ^W^. m.l) 



We shall work within the one elastic constant approximation. Although it is not har d to include more general 



clastic terms this simplification will not affect the qualitative behaviour. The free energy (II. I) describes a first order 
transition from the isotropic to the nematic phase. 

Equation of motion of the nematic order parameter: The equation of motion for the nematic order parameter 
is 1 

(d t +u - V)Q - S(W, Q) = TH (II.2) 



where T is a collective rotational diffusion constant. The first term on the left-hand side of equation ( UL3 ) is the 
material derivative describing the usual time dependence of a quantity advected by a fluid with velocity u. This is 
generalised by a second term 

S(W, Q) = (£D + n)(Q + 1/3) + (Q + 1/3) (^D - fi) 

-2£(Q + I/3)Tr(QW) (II.3) 

where D = (W + W T ) /2 and f2 = (W — W T )/2 are the symmetric part and the anti-symmetric part respectively of 
the velocity gradient tensor W a /3 = dpita- S(W, Q) appears in the equation of motion because the order parameter 
distribution can be both rotated and stretched by flow gradients. £ is a constant which will depend on the molecular 
details of a given liquid crystal. 
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The term on the right-hand side of equation ([1.2) describes the relaxation of the order parameter towards the 
minimum of the free energy. The molecular field H which provides the driving motion is related to the derivative of 
the free energy by 



H 



_8JF_ 
aQ 



(I/3)Tr 



6JF_ 
SQ 



b (Q 2 - (I/3)TrQ 2 ) - cQTrQ 2 + kV 2 Q. 



(II.4) 



Continuity and Navier-Stokes equations: The fluid momentum obeys the continuity 

d t p + d a pu a = 0, (II. 5) 

where p is the fluid density, and the Navier-Stokes equation 

pd t U a + pUpdpU a = d(3T af 3 + dpCT a p + — {dp((S a p - 3d p P )d 1 u 7 + d a up + dpu a ). (II. 6) 

The form of the equation is not dissimilar to that for a simple fluid. However the details of the stress tensor reflect 
the additional complications of liquid crystal hydrodynamics. There is a symmetric contribution 

1 S^F 
+2£(Qa/3 + -8 a p)Q^H ie - dpQ^ 5Q - (11.7) 



and an antisymmetric contribution 
The pressure P is taken to be 



Tap — Qa^/H^p — H al Q 1 p. (II. 



H 



P = pT - -(VQ) . (II.9) 



An earlier development of liquid crystal hydrodynamics in terms of a tensor order parameter was proposed by 
Doi jn|. The Doi theory is based upon a Smoluchowski evolution equation (similar to the Boltzmann equation for 
translational motion) for the orientational distribution function. The main advantage of the approach is the possibility 
of relating the phenomenological coefficients in the equations of motion to microscopic parameters. One omission is 
the lack of gradient terms in the free energy (but see Jl7t ) . Moreover it is necessary to use closure approximations 
to obtain a tractable set of hydrodynamic equations. The Doi and Beris-Edwards equations are very similar: the 
main difference is in the symmetric contribution to the stress tensor. The Doi theory gives a simpler form which is 
incomplete in that it does not obey Onsager reciprocity. (A similar comment applies to all closure relations that we 
have found in the literature.) 

Hydrodynamic equations for the nematic phase were formulated by Ericksen and Leslie These are widely 

used as the Leslie coefficients provide a useful measure of the viscous properties of the liquid crystal fluid. The 
Beris-Edwards equations reduce to those of Ericksen and Leslie in the uniaxial nematic phase when the magnitude of 
the order parameter remains constant. Hence a limitation of the Ericksen-Leslie theory is that it cannot include the 
hydrodynamics of topological defects. For convenience we list below the relationship between the Leslie coefficients 



and the parameters appearing in the equations of motion (II. 2) and (II.6). An outline of their derivation from the 
Beris-Edwards approach is given in Appendix A. 



ai = 


2 9, 


4<z 2 )£ 2 /r 


(11.10) 


a 2 = 


(~|ff(2 + - 


- q 2 )/r 


(11.11) 


a 3 = 


(-^(2 + «)£H 


-q 2 )/r 


(11.12) 


«4 = 


l(i- q ) 2 e/r 


+ 11 


(11.13) 


«5 = 


(^(4-gK 2 + 


±q(2 + q)0/r 


(11.14) 


as = 


(|</(4 -q)e- 


\q{2 + q)0/T 


(11.15) 
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where q is the magnitude of the nematic order parameter and 77 = prj/3. 

A detailed comparison of the theories of liquid crystal hydrodynamics can be found in Beris and Edwards [0 . 



III. A LATTICE BOLTZMANN ALGORITHM FOR LIQUID CRYSTAL HYDRODYNAMICS 



W e n ow d efine a latt ice Boltzmann algorithm which solves the hydrodynamic equations of motion of a liquid crystal 
(11.2), ( |ll.5| ), and (1L6). Lattice Boltzmann algorithms are defined in terms of a set of continuous variables, usefully 
termed partial distribution functions, which move on a lattice in discrete space and time. They were first developed 
as mean-field versions of cellular automata simulations but can also usefully be viewed as a particular finite-difference 
implementation of the continuum equations of motion Q . 

Lattice Boltzmann approaches have been particularly successful in modeling fluids which evolve to minimise a free 
energy [|. It is not proven why this is the case, but one can surmise that the existence of an H-theorem, which 
governs the approach to equilibrium, helps to enhance the stability of the scheme Jl8|,|l9| . 

The simplest lattice Boltzmann algorithm, which describes the Navier-Stokes equations of a simple fluid, is defined 
in terms of a single set of partial distribution functions which sum on each site to give the density. For liquid crystal 
hydrodynamics this must be supplemented by a second set, which are tens or vari ables, and which are related to the 
tensor o rder p arameter Q. A description of the algorithm is given in Section III A and the continuum limit is taken in 
Section III B . A Chapman- Enskog expansion [^o|| showing how the algorithm reproduces the liquid crystal equations 
of motion follows in Section III C . 



A. The lattice Boltzmann algorithm 



We define two distribution functions, the scalars fi(x) and the symmetric traceless tensors Gi(x) on each lattice 
site x. Each fi, Gi is associated with a lattice vector e*j. We choose a nine- velocity model on a square lattice with 
velocity vectors e$ = (±1, 0), (0, ±1), (±1, ±1), (0, 0). Physical variables are defined as moments of the distribution 
function 



(111.I6) 



The distribution functions evolve in a time step At according to 

At 

fi(x + eiAt,t + At) - fi(x, t) = — [C fi (x, t, {fi}) +C fi (x + eiAt, t + At, {/*})] , 
Gi(x + eiAt,t + At) - Gi(x,t) = 

^ [C Gl (x, t, {G 4 }) + C Gl (x + e t At, t + At, {G|})] . 



(III. 17) 



(111.18) 



This represents free streaming with velocity e*; and a collision step which allows the distribution to relax towards 
equilibrium. /* and G* are first or der app roxim ations to fi(x + tiAt, t + At) and Gi(x + e^Ai, t + At) respectively. 
They are obtained from equations ( [II. L7[ ) and ( [11.18 ) but with /* and G* set to fi and G;. Discretizing in this 
way, which is similar to a predictor-corrector scheme, has the advantages that lattice viscosity terms are eliminated 
to second order and that the stability of the scheme is improved. 

The collision operators are taken to have the form of a single relaxation time Boltzmann equation j^] , together with 
a forcing term 



C fi (x,t,{f i })^--(f i (x,t)-ff q (x,t, {fi})) + Pi (x,t, {fi}), 
T f 

C Gl (x, t, {G 4 }) = --(Gi(x, t) - G e t q (x, t, {Gi})) + M z (x, t, {Gi}), 

Tn 



(111.19) 
(111.20) 



The form of the equations of motion and thermodynamic equilibrium follow from the choice of the moments of the 
equilibrium distributions f^ q and G^ 9 and the driving terms pi and M,;. f^ q is constrained by 



^ ft q e-ia = pu a , fi q e ia eip = -<r a /3 + pu a u/3 



(111.21) 
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where the zeroth and first moments are chosen to impose conservation of mass and momentum. The second moment 
of f eq controls the symmetric part of the stress tensor, whereas the moments of pi 

^pi = 0, ^Pieia = df3T af3l ^p i e ia e il 3 = (111.22) 

i i i 

impose the antisymmetric part of the stress tensor. For the equilibrium of the order parameter distribution we choose 

G- 9 = Q, ^ G i qe 'ia = Qu a , ^ G V e iot e iP = Q,U a U0. (III. 23) 

i i i 

This ensures that the order parameter is convected with the flow. Finally the evolution of the order parameter is 
most conveniently modeled by choosing 

^M i = rH(Q) + S(W,Q) = H, 5^M i e to = (5^M i )u a . (IIL24) 

i i i 

which ensures that th e fluid m inimises its free energy at equilibrium. 

Conditions (111.21)— (111.24) can be satisfied as is usual in lattice Boltzmann schemes by writing the equilibrium 
distribution functions and forcing terms as polynomial expansions in the velocity M 

!t q = A s + B s u a e la + C s u 2 + D s u a upe ia eifj + E saf} e ta eip, 
Gl q = J s + K s u a e ia + L s u 2 + N 8 u a upei a eip, 
Pi = T s dpT a pei a , 

M s = R S + S s u a e ia , (111.25) 

where s = e*i 2 £ {0, 1, 2} identifies separate coefficients for different absolute values of the velocities. A suitable choice 
is 

A 2 = (<r xx + 0yy)/16, A x = 2A 2 , A = p- 12A 2 , 

B 2 = p/12, B\ = AB 2 , 

C 2 = -p/16, C\ = -p/8, C a = -3p/4, 

D 2 = p/8, D 1 = p/2 

E 2X 'X = — (T yi/)/l6, E 2 yy — —E 2xX , E 2xy = E 2yX = O X y / Q , 

E\ XX — ^E 2xXl E\yy = 4:E 2yy , 

Jo = Q, 

K 2 = Q/12, K!=4K 2 , 

L 2 = -Q/16, Lj = -Q/8, L = -3Q/4, 

N 2 = Q/8, N a = Q/2 

T 2 = 1/12, T x - 4T 2 , 

R 2 = H/9, Ri = Ro = R 2 

S 2 =H/12, S!=4S 2 , (111.26) 
where any coefficients not listed are zero. 



B. Continuum limit 



We write down the continuum limit of the lattice Boltzmann evolution equations ( III. 17 ) and ( III. 18] ) showing, in 
particular, that the pr edicto r-corrector form of the collision integral eliminates lattice viscosity effects to second order. 
Consider equation ( III. 17 ). Taylor expanding fi(x + eiAt,t + At) gives 



At 2 

fi(x + e,At, t + At)= fi(x, t) + AtDfi(x, t) + —D 2 /^, t) + 0(At 3 ) 



(111.27) 



where D = <9 t + ei a d a . Similarly, expanding the collision term equation([II.19), 
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C fi (x + eiAt, t + At, {fi + AtC fi (x, t, {/,})}) = C fi (x, t, {fi}) 
AtDC fl (x,t,{f t }) + 0(At 2 ) 



(111.28) 



and substituting into equation ([11.17) gives 



At 



Df i (x',t)=Cf i (x,t,{f i }) ~ ^ {D 2 f z (x,t) -DC fi (x,t, {/,})} +0(At 2 ). 



(111.29) 

We see immediately that 

Dfi(x,t) = C fi (x,t, {fi}) + O(At). (IH.30) 
Using equation( III.3C| ) in the expansion ( III. 29 ) it follows that there are no terms of order At in ( [II. 29 ) and 

Dfi(x, t) = C fi (x, t, {fi}) + 0{At 2 ). (111.31) 



A similar expansion of equation ( 111.18 ) leads to 

DGi(x, t) = C Gi (x, t, {Gi}) + 0(At 2 ). 



(111.32) 



In the standard lattice Boltzmann discretization terms of order At appear in equations ( [11.31 ) and ( III. 32 ). These 
are of similar forms to those which arise from the Chapman-Enskog expansion and have been subsumed into the 
viscosity. However this is not gene rally p ossibl e and it is convenient to use the predictor-corrector form for the 
collision term assumed in equations (III. 19) and ([11.20) to eliminate them at this stage. 



C. Chapman-Enskog expansion 



We can now proceed with a Chapman-Enskog expansion, an expansion of the distribution functions about equi- 
librium, which assumes that successive derivatives are of increasingly high order |po| . The aim is to sho w tha t 
equation ([II. 32) reproduces the evolution e quat ion of the l iquid crystal order parameter (EI.2) and equation ( ]III.31| ) 
the continuity and Navier-Stokes equations ( [I.5| ) and ([L6) to second order in derivatives. Writing 



G, = Gf 5 + G\ x > + G\ 

and substituting into ( EII.32| ) using the form for the collision term ( III.2C| ) gives, to zeroth order 

Gf } =Gf + Tg M t . 

Summing over i and using, from equations ( HII.16D and (|III.23| ), 



,(2) 



E G > = Q 



G? 



(111.33) 



(111.34) 



(111.35) 



shows that the zeroth m oment of appears at first order in the Chapman-Enskog expansion. This is as expected 
because, from equation (III. 24), is related to free energy derivatives which will be zero in equilibrium. The 

first moment will also be first order in derivatives. 



It then follows, from substituting equation (III. 33) into equation (III. 32), that the first and second order deviations 
of the distribution function from equilibrium are 



-r g DG? 



TgM, 



(111.36) 



G?) = t 2 D 2 G? - r 2 M„ 



(111.37) 



Using equation ( [11.36 ) in equation ( [II. 33 ), summing over i and using ( III.35| ), ( [II. 23 ), and ( 111.24 ) gives, to first 
order, 



d t Q + d a (Qu a ) =H + 0(d 2 ) 



(111.38) 



G 



The second order term ( III.37| ) gives, after a lengthy calculation, described in Appendix B, a correction 



(111.39) 



This additional term is a feature common to most lattice Boltzmann models of complex fluids, ft is not known whether 
it has a physical orign, but it is very small in all the cases tested so far and has no effect upon the behaviour of the 
fluid. 

A similar expansion for the partial density distribution functions fi gives the continuity and Navier-Stokes equations. 
Writing 



substituting into (111.31) and using the collision operator (111.19) gives 

/i 0) =/r+w 



(111.40) 



(111.41) 



(111.42) 



(111.43) 



Summi ng /j over i and using the constraints on the moments of fi, f^ q and pi, from equations (III. 16), ([11.21) and 
( 111.22 ) respectively 



^d t p + d a pu a + Tfd a ^ P;e iQ ^ = Tfd t 



T f d a 



d t p + d a pu a + T f d a ^2 

i 

d t pu a + dp ^2 JT&ia^d + Tfd t ^2 



The first term in square brackets is second order in derivatives. Therefore 



d t p + d a pu a + T f d a ^2pie ia = Tfd a d t pu a + dp ^ ft q e ia e ll3 + Tfd t y^e^ + 0(d 

\ i / L i i - 



'')■ 



(111.44) 



(111.45) 



We n ow multiply Eq.( [11.40 ) by e ia and sum over i. Using the constraints ( 111.21 ) and ( III.22j ) and the definitions 
(pL16|) 



d t pu a + dp^2 fi 9e iaeif3 + Tfd t ^2p l e ia 

\ i if 



^Pieia + Tfdt 



dtpU a + dp ^ fi qe iaSi0 + Tfd t ),Pi 



+Tfd f3 

So to first order in derivatives 



d t 22 fi 9e ia e iP + d~f ^2 fi 9e iaei0ei~f + T fd 7 22 Pi e ia e ip e i' 



d t pu a + dp ^2 fV^^P + Tfdt ^2pie ia = ^ 

\ i i / i 



Pi e ia + 0{d 2 ). 



(111.46) 



(111.47) 



Placing (III. 47) into the square brackets in equation (111.45) we obtain the continuity equation flll.q ) to second order 
in derivatives 



(d t p + d a pu a ) =0 + O(d 3 ). 



(111.48) 
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Substituting equation ( [11.47 ) into the first square brackets in equatio n (|III.46[) a nd imp osing the constraints on 
the first moment of the pi and the second moment of the f^ q , equations ( p^II. 22[) and ( III.21 ) , gives 



d t {pu a ) + dp{pu a up) = df3a af3 + d/3T a/ 3 



+r f d 



(111.49) 



showing immediately th at the equation of motion ( II. 6 ) is reproduced to Euler level (first order in derivatives). 
From the definitions ( pTI.2C| ) 



i 

y P«eiae^ei 7 = -(dsT$ a 5p 7 + d S T Sl3 5 ai + d S Ts 7 5 al3 ). 



(111.50) 
(111.51) 



Using equations ( III. 50 ) and ( [11.51 ) the viscous terms in the square brackets in equation ( III. 49 ) can be simplified. 
We assume that the fluid is incompressible, ignore terms of third order in the velocities, and furthermore assume that, 
within these second order terms, the stress tensor can be approximated by minus the equilibrium pressure Pq. We 
consider each term in the square brackets in turn: 



1. The first term can be rewritten as 

dt<7 a /3 = —(dpPo)(dtp)5 a /3 = p{d p Po)d 1 u 1 S af 3 



(111.52) 



where the last step follows using the continuity equation (II. 5). 
2. Rewriting 

d t {pu a up) = d t (pu a )uf3 + u a d t (puf3) 



(111.53) 



and replacing time derivatives with space derivatives using the Euler terms in equation (III. 49) one sees that 
this term is zero, given the assumptions listed above. 



3. Using equation (III. 50) 



(111.54) 



4. From equation ( [11.51) the fourth term is of third order in derivatives and can be neglected. 



Replacing the square brackets in the equation ([11.49) with the contributions from 1 and 3 we obtain the incompressible 
Navier-Stokes equation (11.6). 



IV. NUMERICAL RESULTS 



The primary aim of this paper is to describe the details of a numerical algorithm for simulating liquid crystal 
hydrodynamics. Therefore we restrict ourselves here to presenting a few, brief, test cases, aimed at checking the 
approach. Further numerical applications are listed in the summary of the paper and will be presented in detail 
elsewhere. 

In equilibrium with no flow the free energy (II. 1 ) is minimised. For a generic lyotropic liquid crystal we take 
a = (1 — 7/3) and b = c = 7, where 7 = fyLvija. is Doi's excluded volume parameter jl5|,|]]. (L is the molecular 
aspect ratio, <f> the concentration, and V2 and a are 0(1) geometrical prefactors.) At a — by (27c), or 7 = 2.7 for the 
generic lyotropic, there is a first order transition to the nematic phases and as 7 is increased further the nematic order 
parameter q increases. The variation of q with 7 can be calculated analytically. Agreement with simulation results is 
excellent as shown in Figure 1. 
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Imposing a shear on the system in the nematic phase will act to align the director h eld along the flow gradient 



Assuming a steady-state, homogeneous flow and a uniaxial nematic state, it follows from ( [1.2 ) that the angle between 
the direction of flow and the director, 6*, is given by B 

£ CO s26l = (IV.55) 
2 + q 

The simulations reproduce this relation well as shown in Figure 2 for different values of q and £. 



When there is no solution to equation ( IV.55 ) the director tumbles in the flow or may move out of the plane to form 
a log-rolling state Figure 3 gives an example of this type of behavior, showing the director angle as a function 

of time. 

Olmsted and Goldbart Jl(| have argued that shear stress acts to favour the nematic over the isotropic phase. Hence 
application of shear moves the phase boundary, which extends from the first-order equilibrium transition at zero shear 
along a line of first-order transitions which end at a non-equilibrium critical point. Numerical results for this boundary 
are shown in Figure 4. The results are qualitatively similar to those of Jl6|,|l7j who obtained the phase boundary for 
a slightly different model using an interface stability argument. 

On the coexistence line the liquid crystal prefers to phase separate into shear bands jl(||l7|||, coexisting regions 
of different strain rate running parallel to the shear direction. Such shear banding occurs spontaneously in the 
simulations reported here. An example is shown in Figure 5. 



V. SUMMARY AND DISCUSSION 



In this paper we described in detail a lattice Boltzmann algorithm to simulate liquid crystal hydrodynamics. In 
the continuum limit we recover the Bcris-Edwards formulation within which the liquid crystal equations of motion 
are written in terms of a tensor order parameter. The equations are applicable to the isotropic, uniaxial nematic, and 
biaxial nematic phases. Working within the framework of a variable tensor order parameter it is possible to simulate 
the dynamics of topological defects and non-equilibrium phase transitions between different flow regimes. 

Lattice Boltzmann simulations have worked well for complex fluids where a free energy can be used to define 
thermodynamic equilibrium. However previous work has concentrated on self-assembly with much less attention 
being paid to more complex flow properties. The algorithm described here includes coupling between the order 
parameter and the flow. This allows the investigation of non-Newtonian effects such as shear-thinning and shear- 
banding. Examples are given in Section IV. 

There are many directions for further research opened up by the rich physics inherent in liquid crystal hydrodynamics 
and the generality of the Beris-Edwards equations. For example results for liquid crystals under Poiseuille flow show 
that the director configuration can depend on the sample history as well as the viscous coefficients and thermodynamic 
parameters ]2l| ] . The effect of hydrodynamics on phase ordering is being investigated |22] | and it would be interesting 
to study the pathways by which different dynamic states transform into each other. The addition of an electric 
field to the equations of motion will allow problems relevant to liquid crystal displays to be addressed. Numerical 
investigations are proving vital as the complexity of the equations makes analytic progress difficult. 



VI. APPENDIX A 



We outline how the Beris-Edwards equations reduce to those of Ericksen, Leslie, and Parodi in the uniaxial nematic 
phase when the magnitude of the order parameter remains constant. Hence we obtain expressions for the Leslie 



coefficients in terms of the parameters appearing in the equations of motion (II. 2) and flH.6| ) H. 

Taking ft to represent the order-parameter field the Ericksen-Leslie stress tensor and the equation of motion for the 
order parameter are, respectively . 

a a/3 = oiinanpn^npD^p + a A D a p + a^n^n^D^a 

+ aerian^Dpp + a 2 n (3 N a + a 3 n a Np, (VI. 56) 

hi L = iiNp + j2n a D afi (VI.57) 



together with the relations 



71 = "3 - Q-2, (VI. 58) 

72 = «6 — a 5 = a 2 + a 3 . (VI. 59) 
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Th e sec ond of these, known as Parodi's relation, is a result of Onsager reciprocity. (Note that, following the convention 
in (E^ ) , the stress tensor is written so that in the corresponding Navier-Stokes equation one contracts on the second 
index when taking the divergence.) 



The N a are co-rotational derivatives 



N a = d t n a + updpria 



The molecular field h is given by 



5T 
5n„ 



K EL V 2 n„ 



(VI.60) 



(VI.61) 



where the last line assumes the one-elastic constant approximation and £ is a Lagrange multiplier to impose ft 2 = 1. 

To obtain the Ericksen-Leslie-Parodi equations from the tensor formalism uniaxial symmetry is imposed on the 
order parameter 



Q a p = q(n a ni3 - l/3(5 ai a). 



(VI.62) 



where q is the magnitude of the larg est eig envalue. We first obtain an expression for k el in terms of k and show that 
equation (II.4) reduces to the form (VI.61). Using the chain rule 



EL 



5T_ 

tin, 



5 J- dQ 



a/3 



(VI.63) 



-n 5Q afj dn^ 

Substituting H from equation (IL4) into equation ( VI.63| ) , writing Q in uniaxial form and simplifying gives after some 
algebra 



EL 



(VI.64) 



Terms proportional to n M have been omitted as these will only chan ge the magn itude o f the order parameter and the 
Lagrange multipier £ will adjust to prevent this. Hence comparing ( VI.61 ) and ( VI.64 ) 

k el = 2q 2 K . (VI.65) 



Consider now the equation of motion for the order parameter ( VI. 57 ). Solving the Q-evolution equation (|L2|) for 
H, and writing Q in uniaxial form gives 



rff Q/ 3 = q(npN a + n a Np) - q^(D a7 n 7 np + n a n 7 D 7 p) 

2 2 2 

+ g(a - l)£Ai/3 + 2g ^naUpDj^n^ + -q(l - q)£,5 a pD Jl/ n„n 1 . 

Substituting this into equation ( |VI.63| ) yields, after some algebra, 



(VI.66) 



^ = 2q 2 N„ - -q{q + 2)£n a D 



(VI.67) 



where we have again omitted terms proportional to n^. Comparison to equation (VI. 57) gives 

7i = 2<? 2 /r, 



72 



-~q(q + 2)t/T. 



(VI.68) 
(VI.69) 



Finally w e con sider how the stress tens or maps between the two theories. Using equations ( VI.66] ) and ( VI.62| ) the 
symmetric (II. 7) and antisymmetric (II. 8) parts of the Beris-Edwards stress tensor become, respectively, 



rr Q/3 = q 2 {n a N fj - N a n fj ) - q(q + 2)/3£(n a n 7 .D 7/ a - D a7 n 7 np) 
ro-Q/3 : 



(VI.70) 



y (<7 + 2)(n /3 7V Q + n a Np) + — (4 - q)(D al n 1 n /3 + n a n 1 D 10 ) 



- — (<?" 1) D afi - 



., (j + q - q )^n a n/3D lv n v n^ 



+terms in 8 a i3D lv n u n 1 (VI. 71) 

where we have ignored the final, distortion, term in (II. 7). A comparison of flVI.7CD and ( |VI.71| ) to (|yi.56|) gives the 
Leslie coefficients ( |II.10[ ) ( pl7l5| ). (These agree with the expressions given by Beris and Edwards in M], apart for the 
formula for a\. However the formula for ct\ listed in |23| is the same as that calculated here.) 
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VII. APPENDIX B 



We obtain the second order term (III.3E) in t he Cha pman-Enskog expansion for the equati on of m otion of the order 
parameter. Proceeding as in the derivation of ( III.38| ) but including the second order term ( [11.37 ) gives 



d t Q + d a (Qu a ) - H = T g {<9 t 2 Q + 2d a d t {Qu a ) + d a dp(Qu a up) - <9 t H - <9 Q (Hu Q ) j 



(VII. 72) 



where we have used the definitions (ill. 23) and (III. 24) to perform the sums over i. Equation (III. 38) shows that the 
first, half the second and the fourth term in the curly brackets are together of higher order in derivatives and can be 
eliminated. 

We next note that 



d a d t (Qu a ) = d a ( - — (d t p)u a + (d t Q)u a + —d t (pu a ) 
• p p 



(VII. 73) 



The time derivatives can be replaced by sp acial d erivatives by using equations ( III.48| ), (111.38), and (III. 47) re- 
spectively. Substituting back into equation ( VII. 72] ) and ignoring terms in J2iPi e ia ~ dpT a p that contain an extra 
derivative 



OtQ + da(Qu a ) - H = Tg < d, 



Q, 



Q 



d/3(pU{3)Ua - dadf}(Qup)Ua + <9 Q (Hu Q ) 



-Tg <{ d a [ —{dp{pUaU(}) - dpo af3 ) ) + d a dp(Qu a v,p) - <9 Q (H?v 



Rearranging the derivatives this simplifies to 



d t Q + d a (Qu a ) - H = -Tg < da[ —d a Po 



Q 



(VII.74) 



(VII.75) 
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FIG. 1. Equilibrium order parameter q versus ccj). The points are from a simulation and the line is the analytic result. 
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FIG. 2. £ times the cosine of twice the angle between the director and the flow £cos(2#) versus the magnitude of the order 



parameter q. The points are from simulations and the line is the expected value 3g/(2 + g) from Equation (IV. 55) 
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FIG. 3. The components of the director as a function of time for a system changing from a metastable tumbling state to a 
stable log-rolling state. 
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FIG. 4. Phase diagram in the shear stress Yl xy , effective temperature a pfane. (a is the coefficient of the quadratic term in 



the free energy (tl.l).) 
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FIG. 5. Shear bands for a range of strain rates. The bands are formed by the coexistence of isotropic (darker) and nematic 
states. The variation of the strain rate across the system, scaled bylOOF to make it dimensionless, is also shown. 
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